Compare vulnerability results: species vs FE

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

library(terra)
library(oharac)
library(tidyverse)
library(here)
source(here('common_fxns.R'))

1 Summary

Compare vulnerability results between unweighted means (script 2) and functional-vulnerability-weighted means (script 3).

2 Data

Inputs are mean (and sd?) vulnerability maps from scripts 2 and 3.

3 Methods

Loop over several focal stressors. For each stressor, read in the maps of species mean \(\bar x_{spp}\) and fv-weighted mean \(\bar x_{FE}\) (and similar for coefficient of variance). Calculate the following comparisons:

• Percent difference in FE mean relative to species mean: \((\bar x_{FE} - \bar x) / \bar x\)
• Percent difference in FE coefficient of variance relative to species mean: \((\bar {CV}_{FE} - \bar {CV}_{spp}) / \bar {CV}_{spp}\)

squish_rast <- function(r, qtile = .999) {
  r_vals <- values(r); r_vals <- r_vals[!is.na(r_vals)]
  r_zlim <- max(abs(quantile(r_vals, 1 - qtile)), abs(quantile(r_vals, qtile)))
  values(r)[values(r) > r_zlim] <- r_zlim
  values(r)[values(r) < -r_zlim] <- -r_zlim
  
  return(list(r = r, zlim = r_zlim))
}

vuln_spp_dir <- here_anx('_output/vuln_maps/vuln_maps_by_species')
mean_fs_spp <- list.files(vuln_spp_dir, pattern = '_mean.tif', full.names = TRUE)
# sdev_fs_spp <- list.files(vuln_spp_dir, pattern = '_sdev.tif', full.names = TRUE)

vuln_fe_dir <- here_anx('_output/vuln_maps/vuln_maps_by_funct_entity')
mean_fs_fe <- list.files(vuln_fe_dir, pattern = '_mean.tif', full.names = TRUE)
# sdev_fs_fe <- list.files(vuln_fe_dir, pattern = '_sdev.tif', full.names = TRUE)

focal_strs <- c('biomass_removal', 'bycatch', 
                'wildlife_strike',
                'habitat_loss_degradation',
                'nutrient_pollution', 'ocean_acidification', 
                'sst_rise', 'marine_heat_wave',
                'light_pollution')

for(f in focal_strs) { ### f <- focal_strs[3]
  str_name <- str_replace_all(f, '[^a-z]+', ' ') %>%
    str_squish() %>%
    str_to_title()
  
  message('Calculating difference maps for ', tolower(str_name), '...')
  
  mean_f_spp <- mean_fs_spp[str_detect(basename(mean_fs_spp), f)]
  # sdev_f_spp <- sdev_fs_spp[str_detect(basename(sdev_fs_spp), f)]
  mean_rast_spp <- terra::rast(mean_f_spp)
  # sdev_rast_spp <- raster::raster(sdev_f_spp)
  # cv_rast_spp <- sdev_rast_spp / mean_rast_spp
  
  mean_f_fe <- mean_fs_fe[str_detect(basename(mean_fs_fe), f)]
  # sdev_f_fe <- sdev_fs_fe[str_detect(basename(sdev_fs_fe), f)]
  mean_rast_fe <- terra::rast(mean_f_fe)
  # sdev_rast_fe <- raster::raster(sdev_f_fe)
  # cv_rast_fe <- sdev_rast_fe / mean_rast_fe
  
  ### Set up rasters; trim extreme values
  mean_diff_rast <- (mean_rast_fe - mean_rast_spp) / mean_rast_spp
  mean_diff_squished <- squish_rast(mean_diff_rast, qtile = .999)
  # cv_diff_rast <- (cv_rast_fe - cv_rast_spp) / cv_rast_spp
  # cv_r_sq <- squish_rast(cv_diff_rast, qtile = .999)
  
  m_d_mask_0.10 <- m_d_mask_0.25 <- mean_diff_squished$r
  values(m_d_mask_0.10)[values(mean_rast_spp) < .10 & values(mean_rast_fe) < .10] <- NA
  values(m_d_mask_0.25)[values(mean_rast_spp) < .25 & values(mean_rast_fe) < .25] <- NA
  
  message('Plotting difference maps for ', tolower(str_name), '...')
  
  ### Set up plot for diverging palette, and symmetric z limits around zero
  map_cols <- hcl.colors(palette = 'Red-Green', n = 50, rev = TRUE)
  #  [1] "Blue-Red"      "Blue-Red 2"    "Blue-Red 3"    "Red-Green"     "Purple-Green" 
  #  [6] "Purple-Brown"  "Green-Brown"   "Blue-Yellow 2" "Blue-Yellow 3" "Green-Orange" 
  # [11] "Cyan-Magenta"  "Tropic"        "Broc"          "Cork"          "Vik"          
  # [16] "Berlin"        "Lisbon"        "Tofino"   
  
  plot(mean_diff_squished$r, col = map_cols, 
       main = paste0('% Diff in mean vuln: ', str_name),
       zlim = c(-mean_diff_squished$zlim, mean_diff_squished$zlim), 
       legend = TRUE, axes = FALSE)  
  plot(m_d_mask_0.10, col = map_cols, 
       main = paste0('% Diff in mean vuln masked 0.10: ', str_name),
       zlim = c(-mean_diff_squished$zlim, mean_diff_squished$zlim), 
       legend = TRUE, axes = FALSE)  
  plot(m_d_mask_0.25, col = map_cols, 
       main = paste0('% Diff in mean vuln masked 0.25: ', str_name),
       zlim = c(-mean_diff_squished$zlim, mean_diff_squished$zlim), 
       legend = TRUE, axes = FALSE)  
  # plot(cv_r_sq$r, col = map_cols, main = paste0('% Diff in coef of var: ', str_name),
  #      zlim = c(-cv_r_sq$zlim, cv_r_sq$zlim), 
  #      legend = TRUE, axes = FALSE)
}